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Abstract 

We analyze a new Monte Carlo method which uses transition ma- 
trix in the space of energy. This method gives an efficient reweighting 
technique. The associated artificial dynamics is a constrained random 
walk in energy, producing the result that correlation time is propor- 
tional to the specific heat. 

One of the important application of Monte Carlo method [2| is to com- 
pute very high-dimensional integrals which appear in statistical mechanics. 
The method generates a sequence of states Xq, Xi, X2, ■ ■ ■ , with a transi- 
tion probability T{X — * X') = P(X'\X). If we want that the distribution 
of X follows P eq (X), it is sufficient to require that 

P eq (X) T(X -> X') = P eq {X') T(X' -> X). (1) 

This is known as detailed balance condition. 

We consider a simple classical spin model, the Ising model, as an example 
of Monte Carlo dynamics. The spins take only two possible values and live 
on the sites of a lattice, for example on a square lattice. The total energy 
is a sum of interactions between nearest-neighbor sites. In a single-spin-flip 
dynamics, a Monte Carlo move consists of picking a site at random, and 
flipping the spin with probability f| 
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where <jq is the spin value before the flip, and J2i a i is the sum of spins 
of the nearest neighbors. Another popular choice is the Metropolis rate 
min(l, exp(— AH/hsT)) where AH is the energy increase due to flip. 

The local Monte Carlo dynamics has some common features: (1) The 
algorithm is extremely general. It can be applied to any classical model. 
(2) Each move involves 0(1) operations and 0(1) degrees of freedom. (3) 
The dynamics becomes slow near a critical point, characterized by divergent 
time scale, 

T(xL z , T = T C , z&2, (3) 

where L is system linear size. See ref. Q and references therein for recent 
results on z. 

Cluster algorithms || have very different dynamical characteristics. The 
Swendsen-Wang algorithm uses a mapping from Ising model to a type of 
percolation model. Each Monte Carlo step consists of putting a bond with 
probability 

p(ai,(Tj) = 1 - exp(-J(<Ti<7j + l)/k B T) (4) 

between each pair of nearest neighbors, by ignoring the spins and looking 
only at the bonds, we obtain a percolation configurations of bonds ||. A 
new spin configuration is obtained by assigning to each cluster, including 
isolated sites, a random sign +1 or —1 with equal probability. 

Some of the salient features of cluster algorithms: (1) The algorithm is 
applicable to models containing Ising/Potts symmetry. (2) Computational 
complexity is still of O(l) per spin per Monte Carlo step. (3) Much reduced 
critical slowing down. The dynamical critical exponent z sw in r oc L Zsw is 
roughly 0, 0.3, 0.5, 1, in dimensions 1, 2, 3, and 4 or higher, respectively. In 
addition, Li and Sokal Q showed that r > ac for some constant a, and c is 
the specific heat. 

The transition matrix Monte Carlo || is related to the single-spin-flip 
dynamics in the following sense, but it has a totally different dynamics. A 
single-spin-flip Glauber dynamics of the Ising model is described by 



dPja, t) 
di 



{a'} 
N 

^2[-Wi(ai) + Wi{-<Ti)F^P(a,t), (5) 



where N is the total number of spins, and w is given by Eq. (|2|), and Fi 
is a flip operator such that FiP{. . . , <7j, . . .) = P{. . . , — o"j, . . .). Transition 
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matrix Monte Carlo dynamics is defined by 



dP{E,t) 
dt 



J2T(E,E')P(E',t) 



E' 



(6) 



where P(E, t) is the probability of having energy E at time t, and 



The transition matrix T(E, E') has some general properties: (1) The ma- 
trix is banded alone diagonal. (2) The column sum is zero, J2e T{E, E') = 0, 
due to the conservation of total probability. (3) J2 E ,T(E,E')P eq (E') = 0, 
due to existence of equilibrium distribution. (4) The transition rate satisfies 
detailed balance conditions, T(E', E)P eq (E) = T(E, E')P eq (E'). 

The transition matrix Monte Carlo dynamics || has the following inter- 
esting features: (1) It is a constrained random walk in energy space. (2) 
The transition rates are derived from single-spin-flip dynamics. (3) It has a 
fast dynamics, r oc c, and (4) it suggests a different histogram reweighting 
method. 

The artificial dynamics described by Eq. (||) can be implemented on 
computer in at least two different ways, we'll call them algorithm A and B. 

Algorithm A 

1. Do sufficient number of constant energy (microcanonical) Monte Carlo 
steps, so that the final configuration is totally uncorrelated with the 
initial configuration. This step is equivalent to pick a state a at random 
from all states with energy E. 

2. Do one single-spin-flip canonical Monte Carlo move. 

Clearly, this algorithm is not very efficient computationally, due to step 1. 
However, it will be helpful in understanding the dynamics. 

Algorithm B 

• A direct implementation of Eq. (q), i.e., a random walk in energy with 
a transition rate T(E, E'). 



T (E,E') = -±- £ £ T{o,o>), (7) 



H{<j)=E H(a')=E' 
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Then in algorithm B, we need to know T(E, E') explicitly, this can be 
done numerically by Monte Carlo sampling, from 



T(E + AE,E) =w(AE){N(a,AE)) E , AE ^ 0, (8) 

and w(AE) = ±(1 - t&nh(AE/(2k B T)). N(a, AE) is the number of cases 
that energy is changed by AE from E for the N possible single-spin flips. 

Note that computation of (N(a,AE)) E can be done with any sampling 
technique which ensures equal probability for equal energy. We use canonical 
simulations at selected temperatures to compute the microcanonical average 
{N(a, AE)) E so that the total histogram is roughly flat. Alternative sam- 
pling methods are given in ref. |K|. The transition matrix can be formed 
with any temperature. The equilibrium distribution and thus the density of 
states n(E) = P eq (E) exp^/A^T) is obtained by solving 

Y / T(E,E')P eq (E') = 0, (9) 

E< 

or by solving a set of detailed balance conditions. The above scheme is 
similar in spirit to the histogram method of Ferrenberg and Swendsen fll[| , 
and the method has a close connection with, but different from the broad 



histogram of Oliveira et al [12]. In Fig. 1, we show the determination of 
two-dimensional Ising model average energy on a 64 x 64 lattice. The errors 
are very small on a whole temperature range. 

We have more or less a complete understanding of the transition matrix 
Monte Carlo dynamics through exact results in limiting cases. The transi- 
tion matrix T(E, E') can be computed exactly in one-dimensional chain of 
length L (with periodic boundary condition), by some combinatorial con- 
sideration, as 

(k + l)(2k + l) , x , N 

T k , k+1 = { - ^ HI + 7), (10) 

(L-2fc)(L-2fc-l) 
Tk+i,k = 2( ^ L _ (1-7), (11) 

where 7 = tanh(2 J/kgT). The diagonal terms are computed from the 
relation 

Tk-i,k + Tk t k + Tk+i^k = 0, (12) 

and the rest of the elements T kk ' = if |fc — A;'| > 1. The integer k = 
0, 1, 2, . . . , [L/2\ is related to energy by E/J = —L + 4k. While the eigen 
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spectrum at temperature T = can be computed exactly as = — 2(k + 
l)(2k + 1)/(L — 1), the eigenvalues at T > is obtained only numerically. 
The most important feature is that r oc L, given an unusual dynamical 
critical exponent of z = 1 in one dimension. 

The dynamics in any dimensions in the large size limit [13] obeys a linear 
Fokker-Planck equation: 



at' dx v dx +x ^ x > t - 

where t' and x are properly scaled time and energy. 

E - u N 



(13) 



(Nd) 1 / 2 

and t' = bt with 



, u N = E, (14) 



b= lim — 7 — T(E, E')(E' - E) 2 , (15) 

where no is the average energy per spin and d = ksT 2 c is the reduced 
specific heat per spin. The major consequence of this result is that the 
relaxation times are r n = ac'/n, n = 1, 2, 3, • • •, with some constant a. 

The exact results can be interpreted with intuitive pictures. First, we 
consider the result of r oc L in one dimension as T — > 0. At sufficiently 
low temperatures with a correlation length £ comparable with the system 
size L, only the ground state (all spins up or down) and the first excited 
states (with a kink pair) are important. Let's consider the time scale for 
Eq — > E\. A spin with opposite sign is created with probability exp(— AK) 
from Boltzmann weight, where K = J / (ksT), in each of the canonical move. 
Thus 

exp(4if) i 2 
t oc oc — oc L. (16) 

where K is chosen such that there is about one kink pair, so that £ ~ 
exp(2if) ~ L. 

Similarly, the result of r oc c can be obtained by the following argument. 
The transition matrix Monte Carlo is a random walk constrained in the range 
5E, due to the gaussian distribution nature of the equilibrium distribution 
P eq (E). The width of this distribution is related to the specific heat by 
5E 2 = cNksT 2 . Each walk changes E by 0(1). To change E by 5E, we 
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need 8E 2 moves, invoking the theory on random walks. In units of transition 
matrix Monte Carlo steps, 



r « a oc c. 17 

N v ' 

Part of the work presented in this talk is in collaboration with Tay Tien 
Kiat and Robert H. Swendsen. This work is supported in part by a NUS 
Faculty Research Grant PR950601. 
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Figure 1: The average energy per spin of the two-dimensional Ising model 
on a 64 2 lattice by the transition matrix Monte Carlo reweighting method. 
The insert shows the relative error with respect to the exact result (obtained 
numerically based on [|14|]). The canonical simulations are performed at 
25 temperatures, each with 10 6 Monte Carlo steps with a single-spin-flip 
dynamics. 
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